if (Cobb_Douglas)
    tech='_Cobb_Douglas';        
    load ../temp/share_model_age.txt
    share_model_age=reshape(share_model_age, [2 2 T 3]);
else
    tech='';            
end

% reset share parameters
nnkid(1:N)=1;
nkid_0_5(1:N)=1;

if (quo)    
    pbar_age=zeros(N,T);    
else

    insample_baseline=load(['../temp/insample_quo' tech '.txt']);
    pbar_age=load(['../temp/pbar_age' tech '.txt']);
    
    % prices before the price change
    pbar_age=reshape(pbar_age, [N T]);
end

% calculate asset accumulated from t=1 to t=4 (matters only for reduce_W)
asset=zeros(N,1);
for t=1:4
    age(1:N)=t;
    W_m=W_m_base.*exp([(expr_m0+t) (expr_m0+t).^2]*param_W_m);
    W_f=W_f_base.*exp([(expr_f0+t) (expr_f0+t).^2]*param_W_f);

    set_regressors;
    evaluate;
    
    asset=R*asset+W_m.*(5200-l_m)-p.*c-pbar_X*52;
    asset(married==1)=asset(married==1)+W_f(married==1).*(5200-l_f(married==1));
end

if (~quo)
    if (reduce_W)
        name='W';        
        W_m_base=(1-pct_reduction/100)*W_m_base;
        W_f_base=(1-pct_reduction/100)*W_f_base;
        
        % full lifetime income as of t=5 after the wage change
        W=R*asset+(1-pct_reduction/100)*(PDV_W_m2+PDV_W_f2)*5200;
    elseif (reduce_W_temporary)
        name='W_temporary';
        W_m_base=(1-pct_reduction/100)*W_m_base;
        W_f_base=(1-pct_reduction/100)*W_f_base;

        % full lifetime income as of t=5 after the wage change
        W=R*asset+((PDV_W_m2+PDV_W_f2)-(pct_reduction/100)*(PDV_W_m3+PDV_W_f3))*5200;    
    elseif (reduce_W_constant_full_income)
        name='W_constant_full_income';        
        W_m_base=(1-pct_reduction/100)*W_m_base;
        W_f_base=(1-pct_reduction/100)*W_f_base;
    elseif (reduce_P_c)
        name='P_c';
        P_c=(1-pct_reduction/100)*P_c;
    elseif (reduce_p)
        name='p';        
        p=(1-pct_reduction/100)*p;
    end
end

% start from age 5
t=5;
age(1:N)=t;
W_m=W_m_base.*exp([(expr_m0+t) (expr_m0+t).^2]*param_W_m);
W_f=W_f_base.*exp([(expr_f0+t) (expr_f0+t).^2]*param_W_f);

set_regressors;
evaluate;

if (quo)
    dlmwrite(['../temp/insample_quo' tech '.txt'],insample);    
    insample_baseline=insample;
end

%% outcomes
outcome=zeros(11,2);

for m=1:2
    j=0;
    mp=m-1;

    % expenditure
    j=j+1;
    eval(['outcome(j,m)=mean(pbar_X(married==' int2str(mp) '),''omitnan'');']);

    % price
    j=j+1;
    eval(['outcome(j,m)=mean(pbar(married==' int2str(mp) '),''omitnan'');']);
    
    % quantity (5)
    j=j+1;    
    eval(['outcome(j,m)=mean(tau_m(married==' int2str(mp) '),''omitnan'');']);

    j=j+1;
    eval(['outcome(j,m)=mean(tau_f(married==' int2str(mp) '),''omitnan'');']);
    
    j=j+1;
    eval(['outcome(j,m)=mean(g(married==' int2str(mp) '),''omitnan'');']);
    
    j=j+1;
    eval(['outcome(j,m)=mean(Y_c(married==' int2str(mp) '),''omitnan'');']);
    
    j=j+1;
    eval(['outcome(j,m)=mean(X(married==' int2str(mp) '),''omitnan'');']);    

    % so far there are 7 outcomes
    
    % log achievement at age 6 (age 5 investment)
    j=j+1; % 8
    eval(['outcome(j,m)=mean(qd_logPsi(married==' int2str(mp) '),''omitnan'');']);
    
    j=j+1; % consumption unit (9)
    aux=alpha.*qd_logPsi;
    eval(['outcome(j,m)=mean(aux(married==' int2str(mp) '),''omitnan'');']);
    
    max_j=j;
end    

logPsi=zeros(N,1);

% simulate from age 6 to 12 investment
for t=6:T+1
    
    if (quo)
        % save pbar before price change
        pbar_age(:,t-1)=pbar;
        
        % save shares for t-1        
        if (~Cobb_Douglas)
            for m=1:2
                for e=1:2
                    share_model_age(m,e,t-1,:)=share_model(m,e,:);
                end
            end
        end
    end
    
    % logPsi_t
    logPsi=delta2*logPsi+qd_logPsi; % logPsi_{t}=delta2*logPsi_{t-1}+delta1*logX_{t-1}+logtheta

    % update X_t    
    if (t<=T)
        age(1:N)=t;
        nkid_0_5=(age<=5);
        W_m=W_m_base.*exp([(expr_m0+t) (expr_m0+t).^2]*param_W_m);
        W_f=W_f_base.*exp([(expr_f0+t) (expr_f0+t).^2]*param_W_f);
        
        % update regressors for the share parameters
        set_regressors;
        evaluate;
    end
end

if (quo)
    dlmwrite(['../temp/pbar_age' tech '.txt'],reshape(pbar_age,[N*T 1]));    
    
    if (~Cobb_Douglas)
        dlmwrite('../temp/share_model_age.txt', ...
                 reshape(share_model_age,[2*2*T*3 1]));
    end
end

for m=1:2
    
    j=max_j;
    mp=m-1;

    % log achievement at age 13
    j=j+1; % 10
    eval(['outcome(j,m)=mean(logPsi(married==' int2str(mp) '),''omitnan'');']);
    
    j=j+1; % 11
    aux=alpha.*logPsi; 
    eval(['outcome(j,m)=mean(aux(married==' int2str(mp) ...
          '),''omitnan'');']);    
    
end

if (quo)
    csvwrite(['../temp/quo' tech '.csv'],outcome)
else
    
    % compare with quo outcome
    outcome0=csvread(['../temp/quo' tech '.csv']);
    change=outcome-outcome0; % change
    change(1:7,:)=change(1:7,:)./outcome0(1:7,:)*100; % percentage change
    
    change(8,:)=change(8,:)*100; % 100*log change
    change(10,:)=change(10,:)*100; % 100*log change
    
    % consumption units: age 6 skill, age 5 investment
    change(9,:)=(exp( beta^(T-5+1)*delta2^(T-5)*change(9,:) )-1)*100;
    
    % consumption units: age 13 skill, ages 5-12 investment
    change(11,:)=(exp( (1-beta)/(beta^(5-1-T)-1)*change(11,:) )-1)*100;
    
    % save level and change
    csvwrite(['../temp/' name '_' num2str(pct_reduction) tech '.csv'],outcome)
    csvwrite(['../temp/change_' name '_' num2str(pct_reduction) tech '.csv'],change)
end
